{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Privatizing Histograms\n",
    "\n",
    "Sometimes we want to release the counts of individual outcomes in a dataset.\n",
    "When plotted, this makes a histogram.\n",
    "\n",
    "The library currently has two approaches:\n",
    "1. Known category set `make_count_by_categories`\n",
    "2. Unknown category set `make_count_by`\n",
    "\n",
    "The next code block imports just handles boilerplate: imports, data loading, plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "from opendp.meas import *\n",
    "from opendp.mod import enable_features, binary_search_chain, Measurement, Transformation\n",
    "from opendp.trans import *\n",
    "from opendp.typing import *\n",
    "enable_features(\"contrib\")\n",
    "max_influence = 1\n",
    "budget = (1., 1e-8)\n",
    "\n",
    "# public information\n",
    "col_names = [\"age\", \"sex\", \"educ\", \"race\", \"income\", \"married\"]\n",
    "data_path = os.path.join('.', 'data', 'PUMS_california_demographics_1000', 'data.csv')\n",
    "size = 1000\n",
    "\n",
    "with open(data_path) as input_data:\n",
    "    data = input_data.read()\n",
    "\n",
    "def plot_histogram(sensitive_counts, released_counts):\n",
    "    \"\"\"Plot a histogram that compares true data against released data\"\"\"\n",
    "    import matplotlib.pyplot as plt\n",
    "    import matplotlib.ticker as ticker\n",
    "\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_axes([1,1,1,1])\n",
    "    plt.ylim([0,225])\n",
    "    tick_spacing = 1.\n",
    "    ax.xaxis.set_major_locator(ticker.MultipleLocator(tick_spacing))\n",
    "    plt.xlim(0,15)\n",
    "    width = .4\n",
    "\n",
    "    ax.bar(list([x+width for x in range(0, len(sensitive_counts))]), sensitive_counts, width=width, label='True Value')\n",
    "    ax.bar(list([x+2*width for x in range(0, len(released_counts))]), released_counts, width=width, label='DP Value')\n",
    "    ax.legend()\n",
    "    plt.title('Histogram of Education Level')\n",
    "    plt.xlabel('Years of Education')\n",
    "    plt.ylabel('Count')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Private histogram via `make_count_by_categories`\n",
    "\n",
    "This approach is only applicable if the set of potential values that the data may take on is public information.\n",
    "If this information is not available, then use `make_count_by` instead.\n",
    "It typically has greater utility than `make_count_by` until the size of the category set is greater than dataset size.\n",
    "In this data, we know that the category set is public information:\n",
    "strings consisting of the numbers between 1 and 20.\n",
    "\n",
    "The counting aggregator computes a vector of counts in the same order as the input categories.\n",
    "It also includes one extra count at the end of the vector,\n",
    "consisting of the number of elements that were not members of the category set.\n",
    "\n",
    "You'll notice that `make_base_geometric` has an additional argument that explicitly sets the type of the domain, `D`.\n",
    "It defaults to `AllDomain[int]` which works in situations where the mechanism is noising a scalar.\n",
    "However, in this situation, we are noising a vector of scalars,\n",
    "and thus the appropriate domain is `VectorDomain[AllDomain[int]]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Educational level counts:\n",
      " [33, 14, 38, 17, 24, 21, 31, 51, 201, 60, 165, 76, 178, 54, 24, 13, 0, 0, 0]\n",
      "DP Educational level counts:\n",
      " [33, 11, 38, 17, 23, 22, 32, 50, 201, 63, 165, 77, 178, 53, 24, 10, 1, 0, 0]\n",
      "DP estimate for the number of records that were not a member of the category set: 0\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFdCAYAAADBvF6wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqOklEQVR4nO3de5xVdb3/8ddbwFA0QBwJuThgXhCSUfCSplmeFC3vqYwe044dtLKTlRVmP8Oyk5ZaqScN09QkRUHN1FOaeSkVExAQBfMGRxQBR8W7cfn8/ljfwT3DDMzArL1nzbyfj8d+zN7fdfl81h6Yz17ftfb3q4jAzMzM2reNKp2AmZmZrZsLtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgW4cn6QlJ+1U6j0qSdISkFyS9JWmXDdzXfpIWtlVurYj7PUm/KXfctiapWlJI6lrpXKxYXLCt0CTNl/RvjdpOkvT3+tcRMSwi7lvHfjr6H9ELgNMiYrOIeKzxwnTsb6eCXv/4TgXyrM9njQ8FEfHfEfGlHGI1+Pdi1l511D9OZu2KpK4RsaKCKWwDPLGOdUZExDPlSMbMWs9n2NbhlZ6FS9pd0jRJb0haLOmitNoD6efr6ezy45I2kvR9SQskLZF0raSeJfv9QlpWJ+n/NYozXtJkSddJegM4KcV+WNLrkhZJulTSxiX7C0lfkfS0pDcl/UjStpIeSvneWLp+o2NsMldJH5L0FtAFmCXp2fV4/zaRdLWk1yQ9CezWaHlI+mjJ66slnVvy+jBJM9MxPCtpdGr/oqS56Vifk3RKau8B/C+wdcnZ/tbpPb2uZL+Hpssdr0u6T9LQkmXzJZ0habakZZImSeq+Hse+o6S7Jb0q6SlJx6T2PSS9LKlLybpHSJqdnm8kaVw63rr0u9uitfHNSrlgW2fzS+CXEfFhYFvgxtS+b/rZK3UbPwyclB6fAoYAmwGXAkjaCfgVcDzQD+gJ9G8U6zBgMtALmAisBL4BbAl8HNgf+EqjbQ4ERgJ7At8BJgD/DgwEhgO1zRxXk7lGxPsRsVlaZ0REbNvsO9O8H5C9V9um/E5s6YaSdgeuBb5N9j7sC8xPi5cAnwM+DHwR+LmkXSPibeAg4KX0u9gsIl5qtN/tgeuB04Eq4E7gj40+0BwDjAYGAzuTvT8tlj443A38HtgKGAP8StJOEfEI8Dbw6ZJNjkvrAnwNOBz4JLA18BrwP62Jb9aYC7Z1BLems6zXJb1OVkibsxz4qKQtI+KtiJi6lnWPBy6KiOci4i3gTGCMsuvcnwf+GBF/j4h/AWcDjQfmfzgibo2IVRHxbkRMj4ipEbEiIuYDvyb7g17qpxHxRkQ8AcwB7krxl5GddTZ3w9jacm2pGaXvo6QDU/sxwI8j4tWIeAG4uBX7PBm4KiLuTu/DixExDyAi7oiIZyNzP3AXsE8L93sscEfa73Kya/SbAHuVrHNxRLwUEa8CfwRqWpE3ZB8m5kfEb9Pv7DFgCnB0Wn496QOUpM2Bg1MbwKnAWRGxMCLeB8YDn2/l78OsARds6wgOj4he9Q/WPGstdTKwPTBP0qOSPreWdbcGFpS8XkB230fftOyF+gUR8Q5Q12j7F0pfSNpe0u2pK/UN4L/JzrZLLS55/m4TrzejaWvLtaV2LX0fI+LPJfsuPZYFTWzbnIFAk93wkg6SNDV1N79OVvAavx/NaXC8EbEq5Vjay/FyyfN3aP69a842wB6NPgweD3wkLf89cKSkDwFHAjMiYkHJtreUbDeXrIelNb8PswZcsK1TiYinI6KWrIvzfGBy6vpsatq6l8j+8NYbBKwgK6KLgAH1CyRtAvRpHK7R68uAecB2qUv+e4DW/2hanOuGWkRWeEv3XeodYNOS1x8pef4CWVd6A6nITSE7M+6bPmjdyQfvx7qmEWxwvJKUcnxxHdu1xgvA/Y0+xGwWEV8GiIgnyT40HETD7vD6bQ9qtG33iGjL/KyTccG2TkXSv0uqSmdkr6fmVcDS9HNIyerXA9+QNFjSZmRnxJPS3d6TgUMk7ZWum45n3cV3c+AN4C1JOwJfbqPDWleuG+pG4ExJvSUNILs+W2omcJykLumGstJu/iuBL0raP92I1T8d+8bAh8je9xWSDgIOKNluMdBHJTf5NZHTZ9N+uwHfAt4HHlrPY5Sk7qUP4HZge0knSOqWHruV3txGVqS/TnZt/qaS9suBH0vaJu28StJh65mbGeCCbZ3PaOAJZXdO/xIYk64vvwP8GHgwdWPuCVwF/I7sDvLngfdIxSpdY/4acAPZGehbZDdRvb+W2GeQnYm9CVwBTGrD42o211aYpYbfw/5Faj+H7EzyebLrzL9rtN3XgUPIPgAdD9xavyAi/kG6oQxYBtwPbBMRbwL/RVZ4XyN7X24r2W4e2YeQ59LvY+vSgBHxFNnNeJcAr6T4h6T7CdbHXmSXHBo/DiC72ewlsi7288k+aNS7nuwDyl8j4pWS9l+m47lL0pvAVGCP9czNDABFrKvnyczWJZ3Vvk7W3f18hdMxsw7IZ9hm60nSIZI2TdfALwAe54OvLJmZtancCrakgZLulfRkGtzg66n9Z5LmpQENbpHUK7VXS3pX2QALMyVdnlduZm3kMLKu0peA7ci6191lZWa5yK1LXFI/oF9EzEjfUZxONpDAALLrPSsknQ8QEd+VVA3cHhHDc0nIzMyswHI7w46IRRExIz1/k+x7iP0j4q6SO1enUvLVGDMzM2taWa5hp7PnXYBHGi36D7LRm+oNlvSYpPsltXTEIzMzsw4v92Hy0t2zU4DTI+KNkvazyAZ2mJiaFgGDIqJO0kiy4SaHlW6TthsLjAXo0aPHyB133DHvQzAzM2tT06dPfyUiqlqzTa5f60oDGtwO/DkiLippPwk4Bdg/ff+1qW3vA86IiGnN7X/UqFExbVqzi83MzNolSdMjYlRrtsnzLnGRjXI0t1GxHk02C9GhpcU6jQTUJT0fQnbX7XN55WdmZlYkeXaJ7w2cADwuaWZq+x7ZTD8fAu7OajpTI+JUsqH9fihpOdkQkaemWXbMzMw6vdwKdkT8nabHVr6zmfWnkF3rNjMzs0Y8N6uZWSe3fPlyFi5cyHvvvVfpVDqc7t27M2DAALp167bB+3LBNjPr5BYuXMjmm29OdXU16VKltYGIoK6ujoULFzJ48OAN3p/HEjcz6+Tee+89+vTp42LdxiTRp0+fNuu5cME2MzMX65y05fvqgm1mZhVVV1dHTU0NNTU1fOQjH6F///6rX//rX+s7xfkHzjnnHM4888wGbTNnzmTo0KHNbjN+/HguuOCCDY7dlnwN28zMGqged0eb7m/+eZ9d6/I+ffowc+ZMICuUm222GWecccbq5StWrKBr1/UvV7W1tYwePZqf/OQnq9tuuOEGamtr13ufleAzbDMza3dOOukkTj31VPbYYw++853vrHHGO3z4cObPnw/Addddx+67705NTQ2nnHIKK1eubLCv7bffnt69e/PIIx9MZ3HjjTdSW1vLFVdcwW677caIESM46qijeOedNQff3G+//agfVfOVV16huroagJUrV/Ltb3+b3XbbjZ133plf//rXbfwuNOSCbWZm7dLChQt56KGHuOiii5pdZ+7cuUyaNIkHH3yQmTNn0qVLFyZOnLjGerW1tdxwww0ATJ06lS222ILtttuOI488kkcffZRZs2YxdOhQrrzyyhbnd+WVV9KzZ08effRRHn30Ua644gqef/751h9oC7lL3MzM2qWjjz6aLl26rHWde+65h+nTp7PbbrsB8O6777LVVlutsd6xxx7LXnvtxYUXXtigO3zOnDl8//vf5/XXX+ett97iwAMPbHF+d911F7Nnz2by5MkALFu2jKeffrpNvsLVFBdsMzNrl3r06LH6edeuXVm1atXq1/VflYoITjzxxAbXp5sycOBABg8ezP3338+UKVN4+OGHgazr/dZbb2XEiBFcffXV3HfffWtsWxq79CtaEcEll1zSqiK/IdwlbmZm7V51dTUzZswAYMaMGau7nvfff38mT57MkiVLAHj11VdZsGBBk/uora3lG9/4BkOGDGHAgAEAvPnmm/Tr14/ly5c32ZVeH3v69OkAq8+mAQ488EAuu+wyli9fDsA///lP3n777TY42qa5YJuZWbt31FFH8eqrrzJs2DAuvfRStt9+ewB22mknzj33XA444AB23nlnPvOZz7Bo0aIm93H00UfzxBNPNLg7/Ec/+hF77LEHe++9NzvuuGOT251xxhlcdtll7LLLLrzyyiur27/0pS+x0047seuuuzJ8+HBOOeUUVqxY0YZH3VCu82HnzfNhm5ltuLlz5671O8m2YZp6f9vVfNhmZmbWdlywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMwqrkuXLtTU1DBs2DBGjBjBhRdeuHp0sfvuu4+ePXtSU1PD0KFDOeecc9bYfsiQITz11FMN2k4//XTOP//8ZmNWV1c3+F51e+ehSc3MrKHxPdt4f8vWucomm2yyeorNJUuWcNxxx/HGG2+sLs777LMPt99+O2+//TY1NTUccsgh7Lrrrqu3HzNmDDfccAM/+MEPAFi1ahWTJ0/mwQcfbNtjqSCfYZuZWbuy1VZbMWHCBC699FIaD+7Vo0cPRo4cyTPPPNOgvba2lkmTJq1+/cADD7DNNtuwzTbbcPjhhzNy5EiGDRvGhAkT1og3f/58hg8fvvr1BRdcwPjx4wF49tlnGT16NCNHjmSfffZh3rx5bXikreOCbWZm7c6QIUNYuXLl6jHC69XV1TF16lSGDRvWoP1jH/sYG220EbNmzQJoMCPXVVddxfTp05k2bRoXX3wxdXV1Lc5j7NixXHLJJUyfPp0LLriAr3zlKxt4ZOsvty5xSQOBa4G+QAATIuKXkrYAJgHVwHzgmIh4TZKAXwIHA+8AJ0XEjLzyMzOz4vjb3/7GLrvswkYbbcS4cePWKNjwwZzXw4YN49Zbb13dnX7xxRdzyy23APDCCy/w9NNP06dPn3XGfOutt3jooYc4+uijV7e9//77bXRErZfnNewVwLciYoakzYHpku4GTgLuiYjzJI0DxgHfBQ4CtkuPPYDL0k8zM+tknnvuObp06cJWW23F3LlzV1/DXpsxY8ZwwAEH8MlPfpKdd96Zvn37ct999/GXv/yFhx9+mE033ZT99tuvwRSZ0PzUnatWraJXr16rr61XWm5d4hGxqP4MOSLeBOYC/YHDgGvSatcAh6fnhwHXRmYq0EtSv7zyMzOz9mnp0qWceuqpnHbaaWSdry2z7bbbsuWWWzJu3LjV3eHLli2jd+/ebLrppsybN4+pU6eusV3fvn1ZsmQJdXV1vP/++6s/GHz4wx9m8ODB3HTTTUA2/3V9l3sllOUucUnVwC7AI0DfiKif++xlsi5zyIr5CyWbLUxtTc+TZmaFUj3ujlatP7/7ca0L0II7ka39evfdd6mpqWH58uV07dqVE044gW9+85ut3k9tbS3jxo3jyCOPBGD06NFcfvnlDB06lB122IE999xzjW26devG2Wefze67707//v0bTLM5ceJEvvzlL3PuueeyfPlyxowZw4gRI9b/QDdA7tNrStoMuB/4cUTcLOn1iOhVsvy1iOgt6XbgvIj4e2q/B/huRExrtL+xwFiAQYMGjWxuonIza19csNsvT6+Zr0JMrympGzAFmBgRN6fmxfVd3eln/S2ALwIDSzYfkNoaiIgJETEqIkZVVVXll7yZmVk7klvBTnd9XwnMjYiLShbdBpyYnp8I/KGk/QvK7AksK+k6NzMz69TyvIa9N3AC8Likmante8B5wI2STgYWAMekZXeSfaXrGbKvdX0xx9zMzMwKJbeCna5FN3d73/5NrB/AV/PKx8zMmhcRrboj21qmLe8T80hnZmadXPfu3amrq2vT4mJZsa6rq6N79+5tsj9P/mFm1skNGDCAhQsXsnTp0kqn0uF0796dAQMGtMm+XLDNzDq5bt26MXjw4EqnYevgLnEzM7MCcME2MzMrABdsMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MCcME2MzMrABdsMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MCcME2MzMrABdsMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MCcME2MzMrgK557VjSVcDngCURMTy1TQJ2SKv0Al6PiBpJ1cBc4Km0bGpEnJpXbmZmtmGqx93R4nXndz+udTsfv6yV2XQOuRVs4GrgUuDa+oaIOLb+uaQLgdLfyrMRUZNjPmZmZoWVW8GOiAfSmfMaJAk4Bvh0XvHNzMw6kkpdw94HWBwRT5e0DZb0mKT7Je1TobzMzMzapTy7xNemFri+5PUiYFBE1EkaCdwqaVhEvNF4Q0ljgbEAgwYNKkuyZmZmlVb2M2xJXYEjgUn1bRHxfkTUpefTgWeB7ZvaPiImRMSoiBhVVVVVjpTNzMwqrhJd4v8GzIuIhfUNkqokdUnPhwDbAc9VIDczM7N2KbeCLel64GFgB0kLJZ2cFo2hYXc4wL7AbEkzgcnAqRHxal65mZmZFU2ed4nXNtN+UhNtU4ApeeViZmZWdB7pzMzMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzAqgUpN/mJl1SNXj7mjV+vO7H9fylccva2U21pH4DNvMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKILeCLekqSUskzSlpGy/pRUkz0+PgkmVnSnpG0lOSDswrLzMzsyLK8wz7amB0E+0/j4ia9LgTQNJOwBhgWNrmV5K65JibmZlZoeRWsCPiAeDVFq5+GHBDRLwfEc8DzwC755WbmZlZ0VTiGvZpkmanLvPeqa0/8ELJOgtTm5mZmVH+gn0ZsC1QAywCLmztDiSNlTRN0rSlS5e2cXpmZmbtU1kLdkQsjoiVEbEKuIIPur1fBAaWrDogtTW1jwkRMSoiRlVVVeWbsJmZWTtR1oItqV/JyyOA+jvIbwPGSPqQpMHAdsA/ypmbmZlZe9Y1rx1Luh7YD9hS0kLgB8B+kmqAAOYDpwBExBOSbgSeBFYAX42IlXnlZmZmVjS5FeyIqG2i+cq1rP9j4Md55WNmZlZkHunMzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAcivYkq6StETSnJK2n0maJ2m2pFsk9Urt1ZLelTQzPS7PKy8zM7MiyvMM+2pgdKO2u4HhEbEz8E/gzJJlz0ZETXqcmmNeZmZmhZNbwY6IB4BXG7XdFREr0supwIC84puZmXUklbyG/R/A/5a8HizpMUn3S9qnUkmZmZm1R10rEVTSWcAKYGJqWgQMiog6SSOBWyUNi4g3mth2LDAWYNCgQeVK2czMrKLKfoYt6STgc8DxEREAEfF+RNSl59OBZ4Htm9o+IiZExKiIGFVVVVWmrM3MzCqrrAVb0mjgO8ChEfFOSXuVpC7p+RBgO+C5cuZmZmbWnrWoYEvauyVtjZZfDzwM7CBpoaSTgUuBzYG7G319a19gtqSZwGTg1Ih4tan9mpmZdUYtvYZ9CbBrC9pWi4jaJpqvbGbdKcCUFuZiZmbW6ay1YEv6OLAXUCXpmyWLPgx0yTMxMzMz+8C6zrA3BjZL621e0v4G8Pm8kjIzM7OG1lqwI+J+4H5JV0fEgjLlZGZmZo209Br2hyRNAKpLt4mIT+eRlJmZmTXU0oJ9E3A58BtgZX7pmJmZWVNaWrBXRMRluWZiZmZmzWrpwCl/lPQVSf0kbVH/yDUzMzMzW62lZ9gnpp/fLmkLYEjbpmNmZmZNaVHBjojBeSdiZmZmzWtRwZb0habaI+Latk3HzMzMmtLSLvHdSp53B/YHZgAu2GZmZmXQ0i7xr5W+ltQLuCGPhMzMzGxN6zu95tuAr2ubmZmVSUuvYf+R7K5wyCb9GArcmFdSZmZm1lBLr2FfUPJ8BbAgIhbmkI+ZmZk1oUVd4mkSkHlkM3b1Bv6VZ1JmZmbWUIsKtqRjgH8ARwPHAI9I8vSaZmZmZdLSLvGzgN0iYgmApCrgL8DkvBIzMzOzD7T0LvGN6ot1UteKbc3MzGwDtfQM+0+S/gxcn14fC9yZT0pmZmbW2FoLtqSPAn0j4tuSjgQ+kRY9DEzMOzkzMzPLrOsM+xfAmQARcTNwM4Ckj6Vlh+SYm5mZmSXrug7dNyIeb9yY2qpzycjMzMzWsK6C3WstyzZZ184lXSVpiaQ5JW1bSLpb0tPpZ+/ULkkXS3pG0mxJu7boCMzMzDqBdRXsaZL+s3GjpC8B01uw/6uB0Y3axgH3RMR2wD3pNcBBwHbpMRa4rAX7NzMz6xTWdQ37dOAWScfzQYEeBWwMHLGunUfEA5KqGzUfBuyXnl8D3Ad8N7VfGxEBTJXUS1K/iFi07sMwMzPr2NZasCNiMbCXpE8Bw1PzHRHx1w2I2bekCL8M9E3P+wMvlKy3MLW5YJuZWafX0vmw7wXubevgERGSYt1rfkDSWLIucwYNGtTWKZmZmbVLlRitbLGkfgDpZ/0Iai8CA0vWG5DaGoiICRExKiJGVVVV5Z6smZlZe1CJgn0bcGJ6fiLwh5L2L6S7xfcElvn6tZmZWaalQ5OuF0nXk91gtqWkhcAPgPOAGyWdDCwgm/0LsqFODwaeAd4BvphnbmZmZkWSa8GOiNpmFu3fxLoBfDXPfMzMzIrKM26ZmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRVA13IHlLQDMKmkaQhwNtAL+E9gaWr/XkTcWd7szMzM2qeyF+yIeAqoAZDUBXgRuAX4IvDziLig3DmZmXUa43u2cv1l+eRhrVbpLvH9gWcjYkGF8zAzM2vXyn6G3cgY4PqS16dJ+gIwDfhWRLxWmbTMzIqjetwdLV53fvccE7FcVewMW9LGwKHATanpMmBbsu7yRcCFzWw3VtI0SdOWLl3a1CpmZmYdTiW7xA8CZkTEYoCIWBwRKyNiFXAFsHtTG0XEhIgYFRGjqqqqypiumZlZ5VSyYNdS0h0uqV/JsiOAOWXPyMzMrJ2qyDVsST2AzwCnlDT/VFINEMD8RsvMzMw6tYoU7Ih4G+jTqO2ESuRiZmZWBJX+WpeZmZm1gAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQFUej5sM7P2Y3zPVqy7LL88zJrgM2wzM7MCcME2MzMrABdsMzOzAvA1bDPrsKrH3dGq9ed3zykRszbgM2wzM7MCcME2MzMrABdsMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MC8Pewzcys3Wv1d+rP+2xOmVROxQq2pPnAm8BKYEVEjJK0BTAJqAbmA8dExGuVytHMzAqqNRO5QCEmc6l0l/inIqImIkal1+OAeyJiO+Ce9NrMzKzTq3TBbuww4Jr0/Brg8MqlYmZm1n5U8hp2AHdJCuDXETEB6BsRi9Lyl4G+FcvOrINrzTXBjng90KxoKlmwPxERL0raCrhb0rzShRERqZg3IGksMBZg0KBB5cnUrLPrgNcDzYqmYl3iEfFi+rkEuAXYHVgsqR9A+rmkie0mRMSoiBhVVVVVzpTNzMwqpiIFW1IPSZvXPwcOAOYAtwEnptVOBP5QifzMzMzam0p1ifcFbpFUn8PvI+JPkh4FbpR0MrAAOKZC+ZmZmbUrFSnYEfEcMKKJ9jpg//JnZGZm1r61t691mZmZWRNcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKoJLTa9oGatV8xt2Pa93OPT2imVm74jNsMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MCcME2MzMrgA75ta7WfN0JWvmVJ3/dyczMKqBDFmyzDm18z1au7w+ZZh2BC7ZZO9C6QXByTMTM2i1fwzYzMysAF2wzM7MCcME2MzMrgLIXbEkDJd0r6UlJT0j6emofL+lFSTPT4+By52ZmZtZeVeKmsxXAtyJihqTNgemS7k7Lfh4RF1QgJzMzs3at7AU7IhYBi9LzNyXNBfqXOw8zM7Miqeg1bEnVwC7AI6npNEmzJV0lqXflMjMzM2tfKlawJW0GTAFOj4g3gMuAbYEasjPwC5vZbqykaZKmLV26tFzpmpmZVVRFBk6R1I2sWE+MiJsBImJxyfIrgNub2jYiJgATAEaNGhX5Z2uV1OphZs/7bE6ZmFlnk+sw1+uh7AVbkoArgbkRcVFJe790fRvgCGBOuXOzDqA1w3Z6yE4zK5BKnGHvDZwAPC5pZmr7HlArqQYIYD5wSgVyMzMza5cqcZf43wE1sejOcudi1iZ8Vm9mZeDJP8waaf11q5wSMTMr4aFJzczMCsAF28zMrABcsM3MzArABdvMzKwAfNOZtVjugwj4Dmozs2b5DNvMzKwAXLDNzMwKwAXbzMysAHwNu4205vquJ6gwM7PW8hm2mZlZAbhgm5mZFYALtpmZWQH4GnYltGZ2J/D3k83MzGfYZmZmReCCbWZmVgAu2GZmZgXggm1mZlYALthmZmYF4IJtZmZWAC7YZmZmBeCCbWZmVgAu2GZmZgXQ7gq2pNGSnpL0jKRxlc7HzMysPWhXBVtSF+B/gIOAnYBaSTtVNiszM7PKa1cFG9gdeCYinouIfwE3AIdVOCczM7OKa28Fuz/wQsnrhanNzMysU1NEVDqH1SR9HhgdEV9Kr08A9oiI00rWGQuMTS+HA3PKnOaWwCuO2SFidoZjdMyOE88xO1bMHSJi89Zs0N6m13wRGFjyekBqWy0iJgATACRNi4hR5UvPMTtSzM5wjI7ZceI5ZseKKWlaa7dpb13ijwLbSRosaWNgDHBbhXMyMzOruHZ1hh0RKySdBvwZ6AJcFRFPVDgtMzOzimtXBRsgIu4E7mzh6hPyzMUxO3zMznCMjtlx4jlmx4rZ6njt6qYzMzMza1p7u4ZtZmZmTShswS73EKaSrpK0RFJZvkYmaaCkeyU9KekJSV8vQ8zukv4haVaKeU7eMUtid5H0mKTbyxRvvqTHJc1cn7s11zNmL0mTJc2TNFfSx3OOt0M6vvrHG5JOzznmN9K/nTmSrpfUPc94KebXU7wn8jq+pv7/S9pC0t2Snk4/e5ch5tHpOFdJavM7mpuJ+bP0b3a2pFsk9SpDzB+leDMl3SVp6zzjlSz7lqSQtGVbxWsupqTxkl4s+f958Dp3FBGFe5DdkPYsMATYGJgF7JRzzH2BXYE5ZTrGfsCu6fnmwD/LcIwCNkvPuwGPAHuW6Xi/CfweuL1M8eYDW5YjVknMa4AvpecbA73KGLsL8DKwTY4x+gPPA5uk1zcCJ+V8XPVjMWxKdk/OX4CP5hBnjf//wE+Bcen5OOD8MsQcCuwA3AeMKtNxHgB0Tc/PL9Nxfrjk+X8Bl+cZL7UPJLvheUFb/21o5hjHA2e0Zj9FPcMu+xCmEfEA8GqeMRrFWxQRM9LzN4G55DzqW2TeSi+7pUfuNzlIGgB8FvhN3rEqRVJPsv+0VwJExL8i4vUyprA/8GxELMg5TldgE0ldyYroSznHGwo8EhHvRMQK4H7gyLYO0sz//8PIPoSRfh6ed8yImBsRT7VlnBbEvCu9twBTycbHyDvmGyUve9CGf4fW8rf858B32jJWC2K2SlELdqcawlRSNbAL2Rlv3rG6SJoJLAHujojcYwK/IPuPsqoMseoFcJek6Wn0vLwNBpYCv01d/7+R1KMMceuNAa7PM0BEvAhcAPwfsAhYFhF35RmT7Ox6H0l9JG0KHEzDwZfy1DciFqXnLwN9yxS3kv4D+N9yBJL0Y0kvAMcDZ+cc6zDgxYiYlWecJpyWuv6vaskllaIW7E5D0mbAFOD0Rp86cxERKyOihuxT9O6ShucZT9LngCURMT3POE34RETsSjYz3Fcl7ZtzvK5kXWKXRcQuwNtk3ai5S4MQHQrclHOc3mRnnYOBrYEekv49z5gRMZesm/Yu4E/ATGBlnjGbySMoQ29UJUk6C1gBTCxHvIg4KyIGpninrWv99ZU+6H2PnD8UNOEyYFughuwD7oXr2qCoBXudQ5h2BJK6kRXriRFxczljp+7ae4HROYfaGzhU0nyySxuflnRdzjHrzwaJiCXALWSXWfK0EFhY0mMxmayAl8NBwIyIWJxznH8Dno+IpRGxHLgZ2CvnmETElRExMiL2BV4ju9+jHBZL6geQfi4pU9yyk3QS8Dng+PThpJwmAkfluP9tyT5kzkp/hwYAMyR9JMeYRMTidIK0CriCFvwNKmrB7vBDmEoS2fXOuRFxUZliVtXfASppE+AzwLw8Y0bEmRExICKqyX6Pf42IXM/KJPWQtHn9c7KbanK9+z8iXgZekLRDatofeDLPmCVqybk7PPk/YE9Jm6Z/v/uT3XuRK0lbpZ+DyK5f/z7vmMltwInp+YnAH8oUt6wkjSa7ZHVoRLxTppjblbw8jBz/DkXE4xGxVURUp79DC8lu+H05r5iw+kNevSNoyd+gtrwTrpwPsmtV/yS7W/ysMsS7nqzbYjnZL/TknON9gqyLbTZZN99M4OCcY+4MPJZizgHOLvPvdD/KcJc42bcLZqXHE+X495Pi1gDT0vt7K9C7DDF7AHVAzzId4zlkf1znAL8DPlSGmH8j+/AzC9g/pxhr/P8H+gD3AE+T3Z2+RRliHpGevw8sBv5chpjPkN0zVP93qM3u2F5LzCnp39Bs4I9A/zzjNVo+n7a/S7ypY/wd8Hg6xtuAfuvaj0c6MzMzK4CidombmZl1Ki7YZmZmBeCCbWZmVgAu2GZmZgXggm1mZlYALthmOVPm75IOKmk7WtKfKpTPjml2oMckbdtoWeksZjMlXdzE9tVNzXS0gTnVlM5WJOlQlWEWPrMi8de6zMogDfF6E9mY8F3Jvu8+OiKeXY99dY0PJmNYn1zGkc2+dG4Ty+aTzQL1ylq2ryb7vnybDVubRtIaFRG5DUFpVnQ+wzYrg4iYQzYAxHfJxiy+DjhL2fzjj6XJB+rPXv8maUZ67JXa90vttwFPptHa7lA2d/kcScc2jpnOWqfqg3mMe6ez2NOBL0u6t6X5SxqZYs0CvlrSfpKkS0te3y5pv/R8dDqGWZLuSW27S3o4HfNDyubs3hj4IXBsOqs/tnS/6T35azqOe9KIZki6WtLFaT/PSfp8S4/HrIhcsM3K5xzgOLKxvbuTDcO6O/Ap4GdpmNQlwGcim5jkWKC0S3pX4OsRsT3ZGO8vRcSIdKbbVPf6tcB3I2JnshGVfhARdwKXAz+PiE81k+e9JV3i30htvwW+FhEjWnKgkqrIxkc+Km1zdFo0D9gnsglQzgb+O7Ipcs8GJkVETURMarS7S4Br0nFMbPSe9CMbFfBzwHktyc2sqLpWOgGzziIi3pY0CXgLOAY4RNIZaXF3YBDZ/NGXSqohm3Vq+5Jd/CMink/PHwculHQ+Wff030pjKZt/u1dE3J+arqHls3V9qrRLPI0v3yuyOX0hG1LxoKY2LLEn8EB9vhFRPxdwT+CaNFZ0kM25vi4f54M5rn8H/LRk2a2RTZ7wpKTOML2ldWIu2GbltSo9RHb2+VTpQknjycaIHkHWA/ZeyeK3659ExD8l7Uo2pv65ku6JiB/mnHtTVtCwp677Otb/EXBvRByRroXft4Hx3y95rg3cl1m75i5xs8r4M/C1NKsVknZJ7T2BRems8QSgS1MbS9oaeCcirgN+RqOpOiNiGfCapH1S0wnA/ayHyKZafV3SJ1LT8SWL5wM1kjaSNJAPpgicCuwraXDKd4uS46ufCvekkv28CWzeTAoPkc3kVh/7b82sZ9ah+QzbrDJ+BPwCmC1pI+B5suuwvwKmSPoC2XXpt5vZ/mNk171Xkc0A9OUm1jkRuFzSpsBzwBdbmNu9klam57Mj4gtp26skBXBXyboPptyfJJtKcwZARCyVNBa4OR3fErLpWn9K1iX+feCO0pjAOEkzgZ80yudrwG8lfRtY2orjMOtQ/LUuMzOzAnCXuJmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRXA/wfD43fUhEve2gAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# public information\n",
    "categories = list(map(str, range(1, 20)))\n",
    "\n",
    "histogram = (\n",
    "    make_split_dataframe(separator=\",\", col_names=col_names) >>\n",
    "    make_select_column(key=\"educ\", TOA=str) >>\n",
    "    # Compute counts for each of the categories and null\n",
    "    make_count_by_categories(categories=categories)\n",
    ")\n",
    "\n",
    "noisy_histogram = binary_search_chain(\n",
    "    lambda s: histogram >> make_base_geometric(scale=s, D=VectorDomain[AllDomain[int]]),\n",
    "    d_in=max_influence, d_out=budget[0])\n",
    "\n",
    "sensitive_counts = histogram(data)\n",
    "released_counts = noisy_histogram(data)\n",
    "\n",
    "print(\"Educational level counts:\\n\", sensitive_counts[:-1])\n",
    "print(\"DP Educational level counts:\\n\", released_counts[:-1])\n",
    "\n",
    "print(\"DP estimate for the number of records that were not a member of the category set:\", released_counts[-1])\n",
    "\n",
    "plot_histogram(sensitive_counts, released_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Private histogram via `make_count_by`\n",
    "This approach is applicable when the set of categories is unknown or very large.\n",
    "Any categories with a noisy count less than a given threshold will be censored from the final release.\n",
    "The noise scale influences the epsilon parameter of the budget, and the threshold influences the delta parameter in the budget.\n",
    "\n",
    "`ptr` stands for Propose-Test-Release, a framework for censoring queries for which the local sensitivity is greater than some threshold.\n",
    "Any category with a count sufficiently small is censored from the release.\n",
    "\n",
    "It is sometimes referred to as a \"stability histogram\" because it only releases counts for \"stable\" categories that exist in all datasets that are considered \"neighboring\" to your private dataset.\n",
    "\n",
    "I start out by defining a function that finds the tightest noise scale and threshold for which the stability histogram is (d_in, d_out)-close.\n",
    "You may find this useful for your application."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def make_base_ptr_budget(\n",
    "        preprocess: Transformation,\n",
    "        d_in, d_out,\n",
    "        TK: RuntimeTypeDescriptor) -> Measurement:\n",
    "    \"\"\"Make a stability histogram that respects a given d_in, d_out.\n",
    "\n",
    "    :param preprocess: Transformation\n",
    "    :param d_in: Input distance to satisfy\n",
    "    :param d_out: Output distance to satisfy\n",
    "    :param TK: Type of Key (hashable)\n",
    "    \"\"\"\n",
    "    from opendp.mod import binary_search_param\n",
    "    def privatize(s, t=1e8):\n",
    "        return preprocess >> make_base_ptr(scale=s, threshold=t, TK=TK)\n",
    "    \n",
    "    s = binary_search_param(lambda s: privatize(s=s), d_in, d_out)\n",
    "    t = binary_search_param(lambda t: privatize(s=s, t=t), d_in, d_out)\n",
    "\n",
    "    return privatize(s=s, t=t)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I now use the `make_count_by_ptr_budget` constructor to release a private histogram on the education data.\n",
    "\n",
    "The stability mechanism, as currently written, samples from a continuous noise distribution.\n",
    "If you haven't already, please read about [floating-point behavior in the docs](https://docs.opendp.org/en/latest/user/measurement-constructors.html#floating-point).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Educational level counts:\n",
      " [33, 14, 38, 17, 24, 21, 31, 51, 201, 60, 165, 76, 178, 54, 24, 13, 0, 0, 0, 0]\n",
      "DP Educational level counts:\n",
      " {'13': 180, '6': 21, '3': 38, '14': 53, '8': 51, '11': 165, '7': 32, '10': 61, '5': 25, '12': 74, '15': 24, '9': 201, '1': 32}\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAFdCAYAAADBvF6wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqSklEQVR4nO3de5xVdb3/8ddbwFAwQERCLg6YF4RkFLykaZQnRUtNDWX0mJYdtLKjlRVmP8OjnbTUSj1pmGYmKQpqhp7STLRUTEBAFMwbHFEURMW7cfn8/ljfwT3DDMzArL1nzbyfj8d+zN7fdfl81h6Yz17ftfb3q4jAzMzMWrfNKp2AmZmZbZgLtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgW5sn6XFJIyudRyVJOlLS85LekrT7Ju5rpKTFLZVbM+J+X9Kvyx23pUmqkhSSOlY6FysWF2wrNEkLJf1bvbaTJP299nVEDImIaRvYT1v/I3oRcFpEdI2IR+svTMf+dirotY/vViDP2nzW+VAQEf8dEV/JIVadfy9mrVVb/eNk1qpI6hgRqyqYwvbA4xtYZ1hEPF2OZMys+XyGbW1e6Vm4pL0kzZD0hqSXJV2SVrs//Xw9nV1+XNJmkn4gaZGkpZKuk9StZL9fTMuWS/p/9eKMlzRZ0vWS3gBOSrEfkvS6pCWSLpe0ecn+QtLXJD0l6U1J50naQdKDKd+bStevd4wN5irpQ5LeAjoAcyQ9sxHv3xaSrpX0mqQngD3rLQ9JHy15fa2k80teHyFpdjqGZySNSu1fkjQ/Heuzkk5J7V2A/wW2Kznb3y69p9eX7PfwdLnjdUnTJA0uWbZQ0pmS5kpaIWmSpM4bcey7SLpb0quSnpR0TGrfW9JLkjqUrHukpLnp+WaSxqXjXZ5+d1s3N75ZKRdsa29+AfwiIj4M7ADclNoPSD+7p27jh4CT0uNTwCCgK3A5gKRdgV8CxwN9gG5A33qxjgAmA92BicBq4JvANsDHgQOBr9Xb5mBgOLAP8F1gAvDvQH9gKFDTyHE1mGtEvB8RXdM6wyJih0bfmcb9kOy92iHld2JTN5S0F3Ad8B2y9+EAYGFavBT4HPBh4EvAzyTtERFvA4cAL6bfRdeIeLHefncCbgDOAHoBdwJ/rPeB5hhgFDAQ2I3s/Wmy9MHhbuD3wLbAGOCXknaNiIeBt4FPl2xyXFoX4BvA54FPAtsBrwH/05z4ZvW5YFtbcFs6y3pd0utkhbQxK4GPStomIt6KiOnrWfd44JKIeDYi3gLOAsYou879BeCPEfH3iPgXcA5Qf2D+hyLitohYExHvRsTMiJgeEasiYiHwK7I/6KV+EhFvRMTjwDzgrhR/BdlZZ2M3jK0v16aaVfo+Sjo4tR8D/CgiXo2I54FLm7HPk4FrIuLu9D68EBELACLijoh4JjL3AXcB+zdxv8cCd6T9riS7Rr8FsG/JOpdGxIsR8SrwR6C6GXlD9mFiYUT8Jv3OHgWmAKPT8htIH6AkbQUcmtoATgXOjojFEfE+MB74QjN/H2Z1uGBbW/D5iOhe+2Dds9ZSJwM7AQskPSLpc+tZdztgUcnrRWT3ffROy56vXRAR7wDL623/fOkLSTtJmpq6Ut8A/pvsbLvUyyXP323gdVcatr5cm2qP0vcxIv5csu/SY1nUwLaN6Q802A0v6RBJ01N38+tkBa/++9GYOscbEWtSjqW9HC+VPH+Hxt+7xmwP7F3vw+DxwEfS8t8DR0n6EHAUMCsiFpVse2vJdvPJelia8/swq8MF29qViHgqImrIujgvBCanrs+Gpq17kewPb60BwCqyIroE6Fe7QNIWQM/64eq9vgJYAOyYuuS/D2jjj6bJuW6qJWSFt3Tfpd4Btix5/ZGS58+TdaXXkYrcFLIz497pg9adfPB+bGgawTrHK0kpxxc2sF1zPA/cV+9DTNeI+CpARDxB9qHhEOp2h9due0i9bTtHREvmZ+2MC7a1K5L+XVKvdEb2empeAyxLPweVrH4D8E1JAyV1JTsjnpTu9p4MHCZp33TddDwbLr5bAW8Ab0naBfhqCx3WhnLdVDcBZ0nqIakf2fXZUrOB4yR1SDeUlXbzXw18SdKB6UasvunYNwc+RPa+r5J0CHBQyXYvAz1VcpNfAzl9Nu23E/Bt4H3gwY08RknqXPoApgI7STpBUqf02LP05jayIn062bX5m0varwR+JGn7tPNeko7YyNzMABdsa39GAY8ru3P6F8CYdH35HeBHwAOpG3Mf4Brgd2R3kD8HvEcqVuka8zeAG8nOQN8iu4nq/fXEPpPsTOxN4CpgUgseV6O5NsMc1f0e9s9T+7lkZ5LPkV1n/l297U4HDiP7AHQ8cFvtgoj4B+mGMmAFcB+wfUS8CfwnWeF9jex9ub1kuwVkH0KeTb+P7UoDRsSTZDfjXQa8kuIflu4n2Bj7kl1yqP84iOxmsxfJutgvJPugUesGsg8of42IV0raf5GO5y5JbwLTgb03MjczABSxoZ4nM9uQdFb7Oll393MVTsfM2iCfYZttJEmHSdoyXQO/CHiMD76yZGbWonIr2JL6S7pX0hNpcIPTU/tPJS1IAxrcKql7aq+S9K6yARZmS7oyr9zMWsgRZF2lLwI7knWvu8vKzHKRW5e4pD5An4iYlb6jOJNsIIF+ZNd7Vkm6ECAiviepCpgaEUNzScjMzKzAcjvDjoglETErPX+T7HuIfSPirpI7V6dT8tUYMzMza1hZrmGns+fdgYfrLfoy2ehNtQZKelTSfZKaOuKRmZlZm5f7MHnp7tkpwBkR8UZJ+9lkAztMTE1LgAERsVzScLLhJoeUbpO2GwuMBejSpcvwXXbZJe9DMDMza1EzZ858JSJ6NWebXL/WlQY0mAr8OSIuKWk/CTgFODB9/7WhbacBZ0bEjMb2P2LEiJgxo9HFZmZmrZKkmRExojnb5HmXuMhGOZpfr1iPIpuF6PDSYp1GAuqQng8iu+v22bzyMzMzK5I8u8T3A04AHpM0O7V9n2ymnw8Bd2c1nekRcSrZ0H7/JWkl2RCRp6ZZdszMzNq93Ap2RPydhsdWvrOR9aeQXes2MzOzejw3q5lZO7dy5UoWL17Me++9V+lU2pzOnTvTr18/OnXqtMn7csE2M2vnFi9ezFZbbUVVVRXpUqW1gIhg+fLlLF68mIEDB27y/jyWuJlZO/fee+/Rs2dPF+sWJomePXu2WM+FC7aZmblY56Ql31cXbDMzq6jly5dTXV1NdXU1H/nIR+jbt+/a1//618ZOcf6Bc889l7POOqtO2+zZsxk8eHCj24wfP56LLrpok2O3JF/DNjOzOqrG3dGi+1t4wWfXu7xnz57Mnj0byApl165dOfPMM9cuX7VqFR07bny5qqmpYdSoUfz4xz9e23bjjTdSU1Oz0fusBJ9hm5lZq3PSSSdx6qmnsvfee/Pd7353nTPeoUOHsnDhQgCuv/569tprL6qrqznllFNYvXp1nX3ttNNO9OjRg4cf/mA6i5tuuomamhquuuoq9txzT4YNG8bRRx/NO++sO/jmyJEjqR1V85VXXqGqqgqA1atX853vfIc999yT3XbbjV/96lct/C7U5YJtZmat0uLFi3nwwQe55JJLGl1n/vz5TJo0iQceeIDZs2fToUMHJk6cuM56NTU13HjjjQBMnz6drbfemh133JGjjjqKRx55hDlz5jB48GCuvvrqJud39dVX061bNx555BEeeeQRrrrqKp577rnmH2gTuUvczMxapdGjR9OhQ4f1rnPPPfcwc+ZM9txzTwDeffddtt1223XWO/bYY9l33325+OKL63SHz5s3jx/84Ae8/vrrvPXWWxx88MFNzu+uu+5i7ty5TJ48GYAVK1bw1FNPtchXuBrigm1mZq1Sly5d1j7v2LEja9asWfu69qtSEcGJJ55Y5/p0Q/r378/AgQO57777mDJlCg899BCQdb3fdtttDBs2jGuvvZZp06ats21p7NKvaEUEl112WbOK/KZwl7iZmbV6VVVVzJo1C4BZs2at7Xo+8MADmTx5MkuXLgXg1VdfZdGiRQ3uo6amhm9+85sMGjSIfv36AfDmm2/Sp08fVq5c2WBXem3smTNnAqw9mwY4+OCDueKKK1i5ciUA//znP3n77bdb4Ggb5oJtZmat3tFHH82rr77KkCFDuPzyy9lpp50A2HXXXTn//PM56KCD2G233fjMZz7DkiVLGtzH6NGjefzxx+vcHX7eeeex9957s99++7HLLrs0uN2ZZ57JFVdcwe67784rr7yytv0rX/kKu+66K3vssQdDhw7llFNOYdWqVS141HXlOh923jwftpnZpps/f/56v5Nsm6ah97dVzYdtZmZmLccF28zMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzcys4jp06EB1dTVDhgxh2LBhXHzxxWtHF5s2bRrdunWjurqawYMHc+65566z/aBBg3jyySfrtJ1xxhlceOGFjcasqqqq873q1s5Dk5qZWV3ju7Xw/lZscJUttthi7RSbS5cu5bjjjuONN95YW5z3339/pk6dyttvv011dTWHHXYYe+yxx9rtx4wZw4033sgPf/hDANasWcPkyZN54IEHWvZYKshn2GZm1qpsu+22TJgwgcsvv5z6g3t16dKF4cOH8/TTT9dpr6mpYdKkSWtf33///Wy//fZsv/32fP7zn2f48OEMGTKECRMmrBNv4cKFDB06dO3riy66iPHjxwPwzDPPMGrUKIYPH87+++/PggULWvBIm8cF28zMWp1BgwaxevXqtWOE11q+fDnTp09nyJAhddo/9rGPsdlmmzFnzhyAOjNyXXPNNcycOZMZM2Zw6aWXsnz58ibnMXbsWC677DJmzpzJRRddxNe+9rVNPLKNl1uXuKT+wHVAbyCACRHxC0lbA5OAKmAhcExEvCZJwC+AQ4F3gJMiYlZe+ZmZWXH87W9/Y/fdd2ezzTZj3Lhx6xRs+GDO6yFDhnDbbbet7U6/9NJLufXWWwF4/vnneeqpp+jZs+cGY7711ls8+OCDjB49em3b+++/30JH1Hx5XsNeBXw7ImZJ2gqYKelu4CTgnoi4QNI4YBzwPeAQYMf02Bu4Iv00M7N25tlnn6VDhw5su+22zJ8/f+017PUZM2YMBx10EJ/85CfZbbfd6N27N9OmTeMvf/kLDz30EFtuuSUjR46sM0UmND5155o1a+jevfvaa+uVllvBjoglwJL0/E1J84G+wBHAyLTab4FpZAX7COC6yC5YTJfUXVKftB8za2+ae+NTE25ssmJYtmwZp556KqeddhpZ52vT7LDDDmyzzTaMGzeO008/HYAVK1bQo0cPttxySxYsWMD06dPX2a53794sXbqU5cuX07VrV6ZOncqoUaP48Ic/zMCBA7n55psZPXo0EcHcuXMZNmxYix1rc5TlLnFJVcDuwMNA75Ii/BJZlzlkxfz5ks0WpzYXbLM2oGrcHc1af2HnnBKxVundd9+lurqalStX0rFjR0444QS+9a1vNXs/NTU1jBs3jqOOOgqAUaNGceWVVzJ48GB23nln9tlnn3W26dSpE+eccw577bUXffv2rTPN5sSJE/nqV7/K+eefz8qVKxkzZkzFCnbu02tK6grcB/woIm6R9HpEdC9Z/lpE9JA0FbggIv6e2u8BvhcRM+rtbywwFmDAgAHDG5uo3Mxal+YX7OOaF8Bn2BvN02vmqxDTa0rqBEwBJkbELan5ZUl90vI+QO0tgC8A/Us275fa6oiICRExIiJG9OrVK7/kzczMWpHcCna66/tqYH5EXFKy6HbgxPT8ROAPJe1fVGYfYIWvX5uZmWXyvIa9H3AC8Jik2ant+8AFwE2STgYWAcekZXeSfaXrabKvdX0px9zMzMwKJc+7xP8ONHZ734ENrB/A1/PKx8zMGhcRzboj25qmJe8T80hnZmbtXOfOnVm+fHmLFhfLivXy5cvp3LllvvLgyT/MzNq5fv36sXjxYpYtW1bpVNqczp07069fvxbZlwu2mVk716lTJwYOHFjpNGwD3CVuZmZWAC7YZmZmBeCCbWZmVgAu2GZmZgXggm1mZlYALthmZmYF4IJtZmZWAC7YZmZmBeCCbWZmVgAu2GZmZgXggm1mZlYALthmZmYF4IJtZmZWAC7YZmZmBeCCbWZmVgAu2GZmZgXggm1mZlYAHSudgJmZtXHjuzVz/RX55FFwPsM2MzMrABdsMzOzAsitS1zSNcDngKURMTS1TQJ2Tqt0B16PiGpJVcB84Mm0bHpEnJpXbmZmtmmqxt3R5HUXds4xkXYkz2vY1wKXA9fVNkTEsbXPJV0MlF6oeCYiqnPMx8zMrLByK9gRcX86c16HJAHHAJ/OK76ZmVlbUqlr2PsDL0fEUyVtAyU9Kuk+SftXKC8zM7NWqVJf66oBbih5vQQYEBHLJQ0HbpM0JCLeqL+hpLHAWIABAwaUJVkzM7NKK/sZtqSOwFHApNq2iHg/Ipan5zOBZ4CdGto+IiZExIiIGNGrV69ypGxmZlZxlegS/zdgQUQsrm2Q1EtSh/R8ELAj8GwFcjMzM2uVcivYkm4AHgJ2lrRY0slp0RjqdocDHADMlTQbmAycGhGv5pWbmZlZ0eR5l3hNI+0nNdA2BZiSVy5mZmZF55HOzMzMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzczMCqBSY4mbmRnA+G7NWHfFhtexNssF28ysBVWNu6NZ6y/snFMi1ua4S9zMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwKILeCLekaSUslzStpGy/pBUmz0+PQkmVnSXpa0pOSDs4rLzMzsyLK8wz7WmBUA+0/i4jq9LgTQNKuwBhgSNrml5I65JibmZlZoeRWsCPifuDVJq5+BHBjRLwfEc8BTwN75ZWbmZlZ0VTiGvZpkuamLvMeqa0v8HzJOotTm5mZmVH+gn0FsANQDSwBLm7uDiSNlTRD0oxly5a1cHpmZmatU1kLdkS8HBGrI2INcBUfdHu/APQvWbVfamtoHxMiYkREjOjVq1e+CZuZmbUSZS3YkvqUvDwSqL2D/HZgjKQPSRoI7Aj8o5y5mZmZtWYd89qxpBuAkcA2khYDPwRGSqoGAlgInAIQEY9Lugl4AlgFfD0iVueVm5mZWdHkVrAjoqaB5qvXs/6PgB/llY+ZmVmReaQzMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MCcME2MzMrABdsMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MCcME2MzMrABdsMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MCcME2MzMrABdsMzOzAnDBNjMzKwAXbDMzswJwwTYzMysAF2wzM7MCyK1gS7pG0lJJ80rafippgaS5km6V1D21V0l6V9Ls9Lgyr7zMzMyKKM8z7GuBUfXa7gaGRsRuwD+Bs0qWPRMR1elxao55mZmZFU5uBTsi7gderdd2V0SsSi+nA/3yim9mZtaWVPIa9peB/y15PVDSo5Luk7R/pZIyMzNrjTpWIqiks4FVwMTUtAQYEBHLJQ0HbpM0JCLeaGDbscBYgAEDBpQrZTMzs4oq+xm2pJOAzwHHR0QARMT7EbE8PZ8JPAPs1ND2ETEhIkZExIhevXqVKWszM7PKKmvBljQK+C5weES8U9LeS1KH9HwQsCPwbDlzMzMza82aVLAl7deUtnrLbwAeAnaWtFjSycDlwFbA3fW+vnUAMFfSbGAycGpEvNrQfs3MzNqjpl7DvgzYowlta0VETQPNVzey7hRgShNzMTMza3fWW7AlfRzYF+gl6Vsliz4MdMgzMTMzM/vAhs6wNwe6pvW2Kml/A/hCXkmZmZlZXest2BFxH3CfpGsjYlGZcjIzM7N6mnoN+0OSJgBVpdtExKfzSMrMzMzqamrBvhm4Evg1sDq/dMzMzKwhTS3YqyLiilwzMTMzs0Y1deCUP0r6mqQ+kraufeSamZmZma3V1DPsE9PP75S0BTCoZdMxMzOzhjSpYEfEwLwTMTMzs8Y1qWBL+mJD7RFxXcumY2ZmZg1papf4niXPOwMHArMAF2wzM7MyaGqX+DdKX0vqDtyYR0JmZma2ro2dXvNtwNe1zczMyqSp17D/SHZXOGSTfgwGbsorKTMzM6urqdewLyp5vgpYFBGLc8jHzMzMGtCkLvE0CcgCshm7egD/yjMpMzMzq6tJBVvSMcA/gNHAMcDDkjy9ppmZWZk0tUv8bGDPiFgKIKkX8Bdgcl6JmZmZ2Qeaepf4ZrXFOlnejG3NzMxsEzX1DPtPkv4M3JBeHwvcmU9KZmZmVt96C7akjwK9I+I7ko4CPpEWPQRMzDs5MzMzy2zoDPvnwFkAEXELcAuApI+lZYflmJuZmZklG7oO3TsiHqvfmNqqcsnIzMzM1rGhgt19Pcu22NDOJV0jaamkeSVtW0u6W9JT6WeP1C5Jl0p6WtJcSXs06QjMzMzagQ0V7BmS/qN+o6SvADObsP9rgVH12sYB90TEjsA96TXAIcCO6TEWuKIJ+zczM2sXNnQN+wzgVknH80GBHgFsDhy5oZ1HxP2Squo1HwGMTM9/C0wDvpfar4uIAKZL6i6pT0Qs2fBhmJmZtW3rLdgR8TKwr6RPAUNT8x0R8ddNiNm7pAi/BPROz/sCz5estzi1uWCbmVm719T5sO8F7m3p4BERkmLDa35A0liyLnMGDBjQ0imZmZm1SpUYrexlSX0A0s/aEdReAPqXrNcvtdURERMiYkREjOjVq1fuyZqZmbUGlSjYtwMnpucnAn8oaf9iult8H2CFr1+bmZllmjo06UaRdAPZDWbbSFoM/BC4ALhJ0snAIrLZvyAb6vRQ4GngHeBLeeZmZmZWJLkW7IioaWTRgQ2sG8DX88zHzMysqDzjlpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAXQsd0BJOwOTSpoGAecA3YH/AJal9u9HxJ3lzc7MzKx1KnvBjogngWoASR2AF4BbgS8BP4uIi8qdk5mZWWtX6S7xA4FnImJRhfMwMzNr1cp+hl3PGOCGktenSfoiMAP4dkS8Vpm0zMyKo2rcHU1ed+EFn80xE8tTxc6wJW0OHA7cnJquAHYg6y5fAlzcyHZjJc2QNGPZsmUNrWJmZtbmVPIM+xBgVkS8DFD7E0DSVcDUhjaKiAnABIARI0ZEGfI0M2s7xndr5vor8snDmq2S17BrKOkOl9SnZNmRwLyyZ2RmZtZKVeQMW1IX4DPAKSXNP5FUDQSwsN4yMzOzdq0iBTsi3gZ61ms7oRK5mJmZFUGlv9ZlZmZmTeCCbWZmVgAu2GZmZgXggm1mZlYALthmZmYF4IJtZmZWAC7YZmZmBeCCbWZmVgAu2GZmZgXggm1mZlYALthmZmYF4IJtZmZWAC7YZmZmBVCR2brMzFql8d2ase6K/PIwa4ALtpm1WVXj7mjW+gs755SIWQtwl7iZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkB+HvYZmbW6jX7O/UXfDanTCqnYgVb0kLgTWA1sCoiRkjaGpgEVAELgWMi4rVK5WhmZgXVnFHroBAj11W6S/xTEVEdESPS63HAPRGxI3BPem1mZtbutbYu8SOAken5b4FpwPcqlYyZJW3wbMWsaCpZsAO4S1IAv4qICUDviFiSlr8E9K5YdmZtXHOuCXqMbbPKq2TB/kREvCBpW+BuSQtKF0ZEpGJeh6SxwFiAAQMGlCdTMzOzCqvYNeyIeCH9XArcCuwFvCypD0D6ubSB7SZExIiIGNGrV69ypmxmZlYxFSnYkrpI2qr2OXAQMA+4HTgxrXYi8IdK5GdmZtbaVKpLvDdwq6TaHH4fEX+S9Ahwk6STgUXAMRXKz8zMrFWpSMGOiGeBYQ20LwcOLH9GZmZmrVulv4dtZmZmTeCCbWZmVgAu2GZmZgXggm1mZlYALthmZmYF4IJtZmZWAC7YZmZmBeCCbWZmVgAu2GZmZgXggm1mZlYAlZxe0zZR8+YzPq55Ox+/opnZmJlZnnyGbWZmVgAu2GZmZgXggm1mZlYALthmZmYF4IJtZmZWAC7YZmZmBdAmv9bVnK87ASy84LM5ZWJmZtYy2mTBbrbx3Zqxrr+fbBXWnH+v4H+zZm2EC7ZZK9C8QXByTMTMWi1fwzYzMysAF2wzM7MCcME2MzMrgLIXbEn9Jd0r6QlJj0s6PbWPl/SCpNnpcWi5czMzM2utKnHT2Srg2xExS9JWwExJd6dlP4uIiyqQk5mZWatW9oIdEUuAJen5m5LmA33LnYeZmVmRVPQatqQqYHfg4dR0mqS5kq6R1KNymZmZmbUuFfsetqSuwBTgjIh4Q9IVwHlApJ8XA19uYLuxwFiAAQMGlC9hKwYPgmNmbVRFCrakTmTFemJE3AIQES+XLL8KmNrQthExAZgAMGLEiMg/W6ukZg8z60FFzKyFtLZhrstesCUJuBqYHxGXlLT3Sde3AY4E5pU7NzMzs43W3GGDm6kSZ9j7AScAj0mandq+D9RIqibrEl8InFKB3MzMzFqlStwl/ndADSy6s9y5mLUIXzc3szLw5B9m9fi6uZm1Rh6a1MzMrABcsM3MzArABdvMzKwAXLDNzMwKwDedWZO1tkEEzMzaE59hm5mZFYALtpmZWQG4YJuZmRWAr2G3kOZc3/W1XTMzay6fYZuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm5mZFYALtpmZWQG4YJuZmRVAqyvYkkZJelLS05LGVTofMzOz1qBVFWxJHYD/AQ4BdgVqJO1a2azMzMwqr1UVbGAv4OmIeDYi/gXcCBxR4ZzMzMwqrrUV7L7A8yWvF6c2MzOzdk0RUekc1pL0BWBURHwlvT4B2DsiTitZZywwNr0cCswrc5rbAK84ZpuI2R6O0THbTjzHbFsxd46IrZqzQce8MtlILwD9S173S21rRcQEYAKApBkRMaJ86TlmW4rZHo7RMdtOPMdsWzElzWjuNq2tS/wRYEdJAyVtDowBbq9wTmZmZhXXqs6wI2KVpNOAPwMdgGsi4vEKp2VmZlZxrapgA0TEncCdTVx9Qp65OGabj9kejtEx2048x2xbMZsdr1XddGZmZmYNa23XsM3MzKwBhS3Y5R7CVNI1kpZKKsvXyCT1l3SvpCckPS7p9DLE7CzpH5LmpJjn5h2zJHYHSY9KmlqmeAslPSZp9sbcrbmRMbtLmixpgaT5kj6ec7yd0/HVPt6QdEbOMb+Z/u3Mk3SDpM55xksxT0/xHs/r+Br6/y9pa0l3S3oq/exRhpij03GukdTidzQ3EvOn6d/sXEm3SupehpjnpXizJd0labs845Us+7akkLRNS8VrLKak8ZJeKPn/eegGdxQRhXuQ3ZD2DDAI2ByYA+yac8wDgD2AeWU6xj7AHun5VsA/y3CMArqm552Ah4F9ynS83wJ+D0wtU7yFwDbliFUS87fAV9LzzYHuZYzdAXgJ2D7HGH2B54At0uubgJNyPq7asRi2JLsn5y/AR3OIs87/f+AnwLj0fBxwYRliDgZ2BqYBI8p0nAcBHdPzC8t0nB8uef6fwJV5xkvt/clueF7U0n8bGjnG8cCZzdlPUc+wyz6EaUTcD7yaZ4x68ZZExKz0/E1gPjmP+haZt9LLTumR+00OkvoBnwV+nXesSpHUjew/7dUAEfGviHi9jCkcCDwTEYtyjtMR2EJSR7Ii+mLO8QYDD0fEOxGxCrgPOKqlgzTy//8Isg9hpJ+fzztmRMyPiCdbMk4TYt6V3luA6WTjY+Qd842Sl11owb9D6/lb/jPguy0Zqwkxm6WoBbtdDWEqqQrYneyMN+9YHSTNBpYCd0dE7jGBn5P9R1lThli1ArhL0sw0el7eBgLLgN+krv9fS+pShri1xgA35BkgIl4ALgL+D1gCrIiIu/KMSXZ2vb+knpK2BA6l7uBLeeodEUvS85eA3mWKW0lfBv63HIEk/UjS88DxwDk5xzoCeCEi5uQZpwGnpa7/a5pySaWoBbvdkNQVmAKcUe9TZy4iYnVEVJN9it5L0tA840n6HLA0ImbmGacBn4iIPchmhvu6pANyjteRrEvsiojYHXibrBs1d2kQosOBm3OO04PsrHMgsB3QRdK/5xkzIuaTddPeBfwJmA2szjNmI3kEZeiNqiRJZwOrgInliBcRZ0dE/xTvtA2tv7HSB73vk/OHggZcAewAVJN9wL14QxsUtWBvcAjTtkBSJ7JiPTEibiln7NRdey8wKudQ+wGHS1pIdmnj05Kuzzlm7dkgEbEUuJXsMkueFgOLS3osJpMV8HI4BJgVES/nHOffgOciYllErARuAfbNOSYRcXVEDI+IA4DXyO73KIeXJfUBSD+Xlilu2Uk6CfgccHz6cFJOE4Gjc9z/DmQfMuekv0P9gFmSPpJjTCLi5XSCtAa4iib8DSpqwW7zQ5hKEtn1zvkRcUmZYvaqvQNU0hbAZ4AFecaMiLMiol9EVJH9Hv8aEbmelUnqImmr2udkN9Xkevd/RLwEPC9p59R0IPBEnjFL1JBzd3jyf8A+krZM/34PJLv3IleStk0/B5Bdv/593jGT24ET0/MTgT+UKW5ZSRpFdsnq8Ih4p0wxdyx5eQQ5/h2KiMciYtuIqEp/hxaT3fD7Ul4xYe2HvFpH0pS/QS15J1w5H2TXqv5Jdrf42WWIdwNZt8VKsl/oyTnH+wRZF9tcsm6+2cChOcfcDXg0xZwHnFPm3+lIynCXONm3C+akx+Pl+PeT4lYDM9L7exvQowwxuwDLgW5lOsZzyf64zgN+B3yoDDH/RvbhZw5wYE4x1vn/D/QE7gGeIrs7fesyxDwyPX8feBn4cxliPk12z1Dt36EWu2N7PTGnpH9Dc4E/An3zjFdv+UJa/i7xho7xd8Bj6RhvB/psaD8e6czMzKwAitolbmZm1q64YJuZmRWAC7aZmVkBuGCbmZkVgAu2mZlZAbhgm+VMmb9LOqSkbbSkP1Uon13S7ECPStqh3rLSWcxmS7q0ge2rGprpaBNzqi6drUjS4SrDLHxmReKvdZmVQRri9WayMeE7kn3ffVREPLMR++oYH0zGsDG5jCObfen8BpYtJJsF6pX1bF9F9n35Fhu2No2kNSIichuC0qzofIZtVgYRMY9sAIjvkY1ZfD1wtrL5xx9Nkw/Unr3+TdKs9Ng3tY9M7bcDT6TR2u5QNnf5PEnH1o+Zzlqn64N5jHuks9gzgK9Kurep+UsanmLNAb5e0n6SpMtLXk+VNDI9H5WOYY6ke1LbXpIeSsf8oLI5uzcH/gs4Np3VH1u63/Se/DUdxz1pRDMkXSvp0rSfZyV9oanHY1ZELthm5XMucBzZ2N6dyYZh3Qv4FPDTNEzqUuAzkU1McixQ2iW9B3B6ROxENsb7ixExLJ3pNtS9fh3wvYjYjWxEpR9GxJ3AlcDPIuJTjeR5b0mX+DdT22+Ab0TEsKYcqKReZOMjH522GZ0WLQD2j2wClHOA/45sitxzgEkRUR0Rk+rt7jLgt+k4JtZ7T/qQjQr4OeCCpuRmVlQdK52AWXsREW9LmgS8BRwDHCbpzLS4MzCAbP7oyyVVk806tVPJLv4REc+l548BF0u6kKx7+m+lsZTNv909Iu5LTb+l6bN1faq0SzyNL989sjl9IRtS8ZCGNiyxD3B/bb4RUTsXcDfgt2ms6CCbc31DPs4Hc1z/DvhJybLbIps84QlJ7WF6S2vHXLDNymtNeojs7PPJ0oWSxpONET2MrAfsvZLFb9c+iYh/StqDbEz98yXdExH/lXPuDVlF3Z66zhtY/zzg3og4Ml0Ln7aJ8d8vea5N3JdZq+YucbPK+DPwjTSrFZJ2T+3dgCXprPEEoENDG0vaDngnIq4Hfkq9qTojYgXwmqT9U9MJwH1shMimWn1d0idS0/ElixcC1ZI2k9SfD6YInA4cIGlgynfrkuOrnQr3pJL9vAls1UgKD5LN5FYb+2+NrGfWpvkM26wyzgN+DsyVtBnwHNl12F8CUyR9key69NuNbP8xsuvea8hmAPpqA+ucCFwpaUvgWeBLTcztXkmr0/O5EfHFtO01kgK4q2TdB1LuT5BNpTkLICKWSRoL3JKObynZdK0/IesS/wFwR2lMYJyk2cCP6+XzDeA3kr4DLGvGcZi1Kf5al5mZWQG4S9zMzKwAXLDNzMwKwAXbzMysAFywzczMCsAF28zMrABcsM3MzArABdvMzKwAXLDNzMwK4P8DM9dnWtmuRGMAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "from opendp.mod import enable_features\n",
    "enable_features(\"floating-point\")\n",
    "\n",
    "preprocess = (\n",
    "    make_split_dataframe(separator=\",\", col_names=col_names) >>\n",
    "    make_select_column(key=\"educ\", TOA=str) >>\n",
    "    make_count_by(MO=L1Distance[float], TK=str, TV=float)\n",
    ")\n",
    "\n",
    "noisy_histogram = make_base_ptr_budget(\n",
    "    preprocess,\n",
    "    d_in=max_influence, d_out=budget,\n",
    "    TK=str)\n",
    "\n",
    "sensitive_counts = histogram(data)\n",
    "released_counts = noisy_histogram(data)\n",
    "# postprocess to make the results easier to compare\n",
    "postprocessed_counts = {k: round(v) for k, v in released_counts.items()}\n",
    "\n",
    "print(\"Educational level counts:\\n\", sensitive_counts)\n",
    "print(\"DP Educational level counts:\\n\", postprocessed_counts)\n",
    "\n",
    "def as_array(data):\n",
    "    return [data.get(k, 0) for k in categories]\n",
    "\n",
    "plot_histogram(sensitive_counts, as_array(released_counts))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}